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Oscillatory Zoning (OZ) is a phenomenon exhibited by many geologically formed crystals. It is 
characterized by quasi periodic oscillations in the composition of a solid solution, caused by self- 
organization. We present a new model for OZ. The growth mechanism applied includes species dif- 
fusion through the solution bulk, particle adsorption, surface diffusion and subsequently desorption 
or incorporation into the crystal. This mechanism, in particular, can provide the synchronization 
effects necessary to reproduce the layered structure of experimentally obtained crystals, lacking in 
other models. We conduct a linear stability analysis combined with numerical simulations. Our 
results reproduce the experimental findings with respect to the patterns formed and a critical su- 
persaturation necessary for OZ to occur. 

PACS numbers: 47.54.-r, 81.10.AJ, 05.65.+b, 82.40.ck 



I. INTRODUCTION 

Oscillatory zoning (OZ) is a phenomenon describing 
repetitive composition variations of binary solid solutions 
along their core-to-rim profile. Traditionally, it was be- 
lieved to be of rare occurrence and its existence was as- 
cribed to variations of external parameters controlling 
the crystal growth, like temperature or concentration 
fluctuations. However, the development of more sophis- 
ticated observation techniques facilitated the detection 
of this phenomenon in all major classes of minerals and 
a wide range of geological environments 0. In addition 
to naturally obtained samples OZ was experimentally re- 
produced in the absence of external fluctuations. Reeder 
et al. |2] were able to grow calcite crystals exhibiting OZ 
of the Mg dopant and Putnis et al. 0- [3- obtained 
end-member zoning in (Ba,Sr)S04 crystals. 

The experimental setup used by Putnis et al. 0> 
is sketched in Fig. ^ It consists of two reservoirs, one 
filled with an aqueous solution of BaCl2/SrCl2 and the 
other with Na2S04. The two reservoirs are connected 
by a column filled with a silica-gel to inhibit convective 
transport. With the beginning of the experiment the re- 
actants from the reservoirs start to diffuse toward each 
other through the column. As the diffusion fields of Ba 2+ 
and Sr 2+ from one reservoir and S04 2 ~ from the other 
reservoir exceed the nucleation threshold product in the 
vicinity of the column center, nuclei form. The solution is 
then strongly supersaturated with respect to the freshly 
generated crystal seeds and the growth commences in a 
layer of few millimeters in width |fj] . After approximately 
one month the experiment was terminated. The obtained 
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crystals exhibited OZ although no external fluctuations 
were imposed on the system. Thus it has been clearly 
shown that OZ can be also attributed to intrinsic mech- 
anisms resulting in spontaneous structure formation 0. 
The wide range of different crystals concerned suggests a 
certain universality of the underlying mechanism. 
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FIG. 1: Experimental setup in which oscillatory zoned crys- 
tals of (Ba,Sr)S04 were synthesized in Refs. (3,3 0- The 
reactants counterdiffuse in the column and (Ba,Sr)S04 crys- 
tals nucleate. The upper window sketches the structure of the 
nucleation zone and the length scales involved. 

The general principle causing OZ is the autocatalytic 
or inhibiting interaction of the substrate with the end 
member concentration in melt or solution 0. If, for 
example, the crystal is rich in component A, this will 
lead to increased growth of this component in an auto- 
catalytic way. Its supply will eventually be limited by 
diffusion. During this phase the disfavored component 
B will accumulate in the solution, leading to a slight in- 
crease of B deposition. However, any small increase in B 
will show autocatalytic effects, whereas the growth of A 
slowly decreases. The combination of a relatively high B 
concentration in combination with the autocatalytically 
increasing growth rate will then lead to a phase of B dom- 
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inated growth. With this, A and B have switched roles 
and one half cycle is completed. If the interaction of those 
two processes is interrupted, for example by stirring no 
OZ will be observed 0. 

The specific interaction of autocatalytic growth and 
component accumulation is subject to the scenario em- 
ployed and gives rise to different nonlinear schemes. The 
first quantitative model derived by Haase et al. 9] de- 
scribes self-organized oscillatory zoning from the melt 
applying moving boundaries and a generically autocat- 
alytic growth term. In a subsequent series of publica- 
tions Wang and Merino introduced the boundary layer 
approximation for the treatment of zoned crystals grown 
hydrothermally |lOj . from solution and from melt 
Hal . Later, non generic growth terms derived from the 
physics of growth processes from the melt were intro- 
duced by L'Heureux u sing constitutional undercool- 
ing and by Wang and Wu [Tj| employing the excess en- 
thalpy of crystallization. 

The most sophisticated models currently available for 
end member OZ from solution have been developed by 
L'Heureux et al. [TH fl6l These models apply the 
boundary layer approximation, as well, and in addition to 
the otherwise deterministic nature the influence of noise 
on the onset of OZ is studied. The non-linear growth 
term applied is phcnomenologically obtained from the lo- 
cal probability to find a matching kink site as proposed 
by Markov ■ However, the local nature of this mech- 
anism does not provide the synchronization effects nec- 
essary to describe homogeneous growth fronts resulting 
in the rin g lik e composition oscillations found in the ex- 
periment \u\ . 

In this paper we present a boundary reaction diffusion 
model for OZ in binary solid solutions grown from aque- 
ous solution. We abandon the boundary layer approxi- 
mation and explicitly treat the diffusion above the crys- 
tal without further approximations. Furthermore, the 
growth rate, being the central ingredient of every model, 
is derived directly from considerations of the physical 
growth mechanisms. We apply the concept of layer-by- 
layer growth under continuous generation of new steps 
which is, e.g., relevant for growth by screw dislocations 
or 2D nucleation. The growth mechanism results as an 
interplay of different processes including bulk diffusion, 
adsorption, surface diffusion, and eventually desorption 
or incorporation into the crystal. The non-linearity nec- 
essary to generate OZ is obtained by the composition- 
dependence of the mean life time of adatoms in the ad- 
sorbed layer or, equivalently, the interaction of adatoms 
with the crystal surface. 

In Sect. II the different aspects of the model are intro- 
duced. The resulting model equations are summarized in 
Sect. III. Then, Sect. IV analyses this model close to the 
stationary point. Sect. V presents its numerical analy- 
sis, including all nonlinear effects. Finally, the obtained 
results are discussed in Sect. VI and concluded in Sect. 
VII. 



II. PHYSICAL BACKGROUND 

The model under consideration describes the diffusion 
processes in the bulk solution, the growth process follow- 
ing from the coupling between crystal and solution, and 
the evolution of the crystal composition. 



A. General 

Based on the slow crystal growth observed in the 
experiments LH 0- we apply screw dislocations as 
the step generating mechanism and describe the subse- 
quent growth process via step advance. Screw dislocation 
driven growth can cross over to two dimensional nucle- 
ation, as shown by Pina et al. 4}. However, this will not 
affect the validity of the model because the specific pro- 
cess of step generation is not of immediate importance. 
We assume that after nucleation the crystal surface will 
reach a steady-state when the density of the step generat- 
ing islands or spirals does not change any more, because 
of coalesence of the terraces according to 19]. A refined 
model would be necessary to take into account any effects 
related to anisotropic growth as found by Pina et al. [2(j 
or to account for the curvature of small spirals. Since we 
are interested in the basic mechanisms we refrain from 
such detailed modelling and just consider infinite step 
trains which on average are a distance I apart; see Fig. [21 
Typically, I is in the nanometer-regime. The coordinates 
z, characterizing the distance from the crystal surface 
and x, orthogonal to the steps, are indicated. This de- 
scribes a one-dimensional crystal surface which starts at 
z = 0, i.e. the total system is a 2D-system. 

In general, quantities like the solute concentration 
C(x, z) depend on x and z. Conceptually, the in- 
dependence can be separated into two different contri- 
butions. First, there exists a periodic contribution with 
quasiperiod /. It expresses the fact that the concentra- 
tion close to the steps will be smaller. However, in the 
limit, considered below, only the concentration, averaged 
over the length scale of I, will enter. Therefore this in- 
dependence of C(x, z) is not relevant. Second, there can 
be variations on length scales much larger than /. Experi- 
mentally, it is observed that the growth behavior does not 
change along the surface of one crystallite of size 150 \x at 
a given time. Furthermore, a straightforward extension 
of the stability analysis, presented below, shows that the 
maximum instability for fluctuations along one crystal 
surface are for zero wave vector. Thus it is realistic to 
hope that the leading mechanism of OZ can be derived 
from study of the z-dependence alone. In the present 
work, we will therefore neglect a possible long-range in- 
dependence and restrict ourselves to the ID model. 

In the present paper we deal with a two species model 
(i = 1,2) for the crystal growth from solution (Ba and 
Sr, respectively), thereby neglecting possible variations 
of the SO,i 2 ~ concentration. This can be justified in two 
ways. First, using the OZ model by L'Heureux [l5T | we 
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FIG. 2: Scheme of the 2D-model. The concentration is av- 
eraged over a sufficiently large x-region, giving rise to a ID- 
model with a concentration C(z). 



have verified that a system with artificially fixed SC>4 2 ~ 
concentration exhibits basically the same dynamics |2l| . 
Second, it describes also the solid state formation for the 
three component system SO4 2- , Ba, and Sr, studied by 
Putnis et al. when the concentration of SC>4 2 ~ is suffi- 
ciently high. The mole fraction of component 1 on the 
crystal surface is denoted by x (0 < x < 1). Through- 
out this work we assume the same molar volume for both 
species in the solid phase. 



B. Crystal-Solution Interface 

The central element of this model is the coupling term 
between the crystal growth rate and the surface species 
concentration. We have used the classical approach by 
Gilmer et al. which itself builds on the BCF-theory 
[l9| : solute particles adsorb on the crystal surface and 
diffuse along it. If before desorbing they come into con- 
tact with a step on the crystal surface, they are incorpo- 
rated. If they do not meet a step in time, they desorb 
and become part of the solution again. 

On their path along the crystal surface the adatoms 
experience many different local environments depending 
on the crystal composition. The mean energy difference 
Ei(x) of the adsorption reflects the interaction with the 
crystal surface and the solvation process. Within the 
mean-field description it can be approximated by a linear 
combination of the different adsorption energies: 



E 1 (x)=xEu + (l-x) E12 
E 2 (x) = xEi2 + (l-x)E 2 2 



sol 



Eg*. 



(1) 



Here, Eij = Eji is the adsorption energy for component i 
on a surface formed by component j and £f ol represents 
the solvation energy of species i. 

In the subsequent mathematical treatment it is useful 
to rewrite expressions Q in the symmetrized form 

Ei(x) = 2 + (1 - X)0 + V) k B T + AE^ , 
E 2 ( X ) =2( X 6+(1- X )4> ~ 1) k B T + AE£ , (2) 



where the dimensionless potentials 

(j>=(E 22 -E n )/4k B T, 
9={2E 12 -(E n +E 22 ))/4k B T, 

ri = (£2°' - £i o1 ) I^bT . 
and the mean homogeneous adsorption energy 
AE^ = (l/2)(E n + E 22 ) - {l/2)(Ef + Ef) 



(3) 



(4) 



have been introduced. Here, k B is the Boltzmann con- 
stant and T the temperature. The potential cf> represents 
the asymmetry between homogeneous adsorption ener- 
gies, whereas 9 is a measure for the preference of homo- 
geneous over heterogeneous adsorption. We consider the 
case < 4> < 0. The potential 4> can be assumed to be 
nonnegative due to symmetry reasons; see Eq. @. The 
limit cf> — 8 implies Ei 2 = E 22 , so that the crystal growth 
properties of species 2 are independent of the composi- 
tion of the crystal surface. The case <f> > 9, corresponding 
to a different type of crystal growth instability, is beyond 
the scope of the present paper. The last parameter 77 
reflects the solution energy difference of the two types of 
particles. 

For the crystal growth two time scales are of primary 
importance. First, the inverse of the adsorption time r a 
denotes the rate with which a particle in the solution 
layer above the crystal surface adsorbs. Thus, aCf/r a 
is the particle flux on the crystal surface where a is the 
typical distance between atoms and C? — Ci(z — 0) the 
mean concentration of component i in the solution just 
above the crystal surface. We assume that r a is the same 
for both species. Second, r^j is the mean residence time 
of adatoms on the surface. Detailed balance between 
both time scales requires 



TdAx) = r a exp[-Ei(x)/k B T] 



(5) 



Thus the composition-dependence of Td,% is due to the 
composition-dependence of the adsorption potentials 
Ei(x)- From Td,i the mean diffusion length /| can be 
obtained using the Einstein relation 



r a exp 



E l ( X )/k B T 



(6) 



The adatom diffusion coefficient D s characterizes the ele- 
mentary atomic movements of the adatoms on the crystal 
surface. We also assume that D s is the same for both the 
species. 

By substitution of equation one explicitly obtains 
for the mean diffusion length 



l\{x)=l D f n exp[-(-x0+(l-x)0)] 

iKx) = IdU 1 ex P [- ( X e + (1 - xM 



with 



l D := \/DsTa exp(-A^ 1 d /2fc B T) 



(7) 



(8) 
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and 



/„ := exp(-r)). 



(9) 



Now the partial growth rate of species i can be ex- 
pressed as a combination of the adsorption flux and a 
success factor qi, describing the fraction of adatoms that 
will actually contribute to crystal growth, whereas the 
others dcsorb: 



■c? 



flux to surface 



Qiix) ■ 



(10) 



Following Ref. pJj we have neglected a term, representing 
the equilibrium concentration which is irrelevant under 
significant supersaturation as present in the experiment. 
The adsorption of a particle is hindered by the breakup 
of the solution shell and consequently takes longer than 
a normal diffusion step in the solution bulk. It is reason- 
able to assume that the adsorption time scale r a is much 
larger than the typical time scale a 2 /D of free diffusion 
in the solution, i.e. y/Dr a ^> a. Since I is also a micro- 
scopic length scale one may even expect that y/D r a 3> I. 
Then the location of successive adsorption processes of 
the same particle are uncorrelated with respect to the po- 
sition of the steps and this simple probabilistic approach 
becomes justified. 

The growth is supported by surface diffusion from both 
the sides of a step. In the limit if <C / the atoms ad- 
sorbing within the distance i| of a step will typically 
contribute to crystal growth, whereas all the remain- 
ing adatoms will desorb after some time. Therefore, the 
probability to meet a step is of the order of If /I and can 
more precisely be expressed as 



Qiix) 



2?f(x) 
I 



(11) 



For arbitrary ratio I* /I one obtains an additional factor 
tanhfZ/2 1*) which in the limit If <C I turns out to be 
unity |22J. 

For later convenience we gather the system constants 
into a characteristic time 



2 l D 



and defining 



nix) 



ID 



the individual growth rates can be written as 



nix,c?) 



nix) 



(12) 



(13) 



(14) 



Then the resulting continuity relation at the crystal sur- 
face reads 



D, 



dCi(z,t) 



dz 



nix,Ci). 



(15) 



z=0 



This expresses the particle flux, on the left side, from 
the solution to the crystal and, on the right side, in the 
continuous solution close to the surface. An analogous 
boundary condition has been used in the classical work 
by Gilmer et al. 



C. Bulk Solution 

The description of the bulk solution follows the follow- 
ing picture: At the beginning of the experiment the com- 
ponents begin to diffuse from the reservoirs into the gel 
column. As their diffusion profiles begin to overlap close 
to the middle of the column and the nucleation barrier is 
overcome, crystallites form. Due to the the narrow nucle- 
ation zone |6( these nuclei must act as an effective sink 
with respect to the current of Ba 2+ , Sr 2+ , and SC>4 2 ~ 
from their reservoirs. 

Based on these considerations, this model is set up as 
a source-sink-system with a gradient in between. In this 
aspect it differs distinctly from the models proposed by 
L'Heureux et al. [Tj, [T(| H3 , where the crystal is consid- 
ered to be growing through a homogeneously supersatu- 
rated medium. Consequentially an analogous boundary 
layer approximation cannot be applied to the present sys- 
tem. 

Experimentally, the following scales are observed 
(i) Growth velocity V ~ 10~ 8 cm/s as estimated from 
the crystal thickness and the total growth time, (ii) 
Thickness H oz ss 10 -5 cm of individual layers, (iii) Time 
T oz « H oz /V w 10 5 s during which one pattern layer is 
formed, (iv) Bulk diffusion constant D rj 10~ 5 cm 2 /s. 

^From these observables two important length scales 
can be estimated, (i) The length scale L oz characterizes 
the spatial variations of the species distribution caused by 
oscillatory zoning. It is given by L oz w (DT oz ) - 5 w 1 cm 
which is much larger than the crystal size. This is con- 
sistent with the fact that the growth behavior does not 
change along one crystal surface, because during times 
of small change in \ the information about the local sur- 
face concentration can spread over the whole region. Be- 
sides, it rationalizes the mean field approach discussed 
below, (ii) Furthermore, one may wonder to which de- 
gree the motion of the solid phase boundary caused by 
crystal growth can affect the diffusion fields. For small 
length scales the concentration field is determined by 
diffusion (w (Dt) - 5 ) whereas for large scales it is de- 
termined by the (nearly) constant growth of the phase 
boundary (« vt). The crossover length scale L v for which 
both processes are equally relevant can be thus estimated 
as L v sa D/V which for the present situation is close to 
10 3 cm. This scale exceeds the system size by orders of 
magnitude. 

Because of this estimate the effect of the growth in- 
duced motion of the crystallites boundaries on the solute 
diffusion is ignorable and the crystal surface can be re- 
garded as fixed with respect to the mathematical model. 

Finally, we choose the external boundary condition 
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such that the reservoir is characterized by a constant in- 
flux Gi of solute into the system at z = L where L is an 
arbitrary large length scale. 



D. Surface composition evolution 

To complete the model, a governing equation for the 
evolution of the crystal surface composition is necessary. 
Assuming a homogeneous distribution of the components 
throughout the surface, the composition change with the 
next time increment dt can be expressed as a function of 
the current composition X an< i the relation of the growth 
rates tv: 



dx 
dt 



(16) 



This relation reflects mass conservation and has already 
been used in previous work on OZ |15|| . 



III. THE BOUNDARY-REACTION-DIFFUSION 
MODEL 

The ID formulation of the model discussed above can 
be described as follows. Diffusion of the components i — 
1, 2 through the solution is considered within the region 
z G [0, L] and is described by the equation 



dd(z,t) 
dt 



= D. 



d 2 C l jz,t) 
dz 2 



(17) 



where L should be chosen large enough to satisfy 
dC/dt = 0. At the external boundary z = L the influx 
of both the components is fixed 



G, = D, 



dCi(z,t) 



dz 



(18) 



At the crystal surface (z — 0) the diffusion flux and the 
rate of the crystal growth are related by (using expres- 
sions and 



D, 



dCi(z,t) 



dz 



nix) 



(19) 



The time scales Tj, as defined in Eq. l|13|l . are given by 



n(x)=Tf- 1 exp{- X (ct> + 9)+6} , 

rtx) = r fv ex P \x(Q - 4>) + 4>] ■ 



(20) 



The constants r and have been defined in Eqs. Ijl2|l 
and ©, respectively. Finally, for the solid composition 
evolution one obtains (compare Eq. (|1<3|) ^ 



dx 
dt 



(i-x) 



aC 9 s 



aCj 

n{x) ^ t 2 (x) 



(21) 



In order to complete the description of OZ it is nec- 
essary to calculate the resulting structure of the crystal, 



1*6. Xcrystal (^crystal)- For this purpose we introduce ((t) 
as the location of the crystal surface at time t in a co- 
ordinate system which does not move with the crystal 
surface. Defining t z as the time for which £(t) = z CT ystai 
one has 



Xcrystal (^crystal) — X(^z) • 

The function £(t) can be easily obtained from 

aC{ aCi " 



dQ/dt 



nix) rzix) 



(22) 



(23) 



where the left side can be interpreted as the time- 
dependent growth velocity proportional to the cumula- 
tive species flux at the surface. 



IV. STABILITY ANALYSIS 

The following two sections analyze the behavior of 
the crystal growth within the boundary-reaction diffu- 
sion model specified in the previous one. Particularly we 
study the crystal growth properties close to the station- 
ary point by means of linear stability analysis. 



A. Stationary solution 

Both, Eqs. (|17|) and (|21|) describe the time-dependence 
of the underlying system. For each equation and for fixed 
X the disappearance of the time-derivative yields some 
ratio Cf/Cf, respectively. 

Evidently, (d/dt)xit) = implies (using /^e* = 1 for 
simplicity) 



X 



■exp(6(l-2 X )) = N x (x). 



C]_ x nix) 

Cf l-xr 2 (x) 

(24) 

This nullcline N x ix) is shown in Fig. |3J It possesses a 
decreasing branch located in the region x € (x-?X+) 



X± := 



1± 



9-2 



(25) 



when the parameter 9 exceeds the critical value 9 C = 2. 
For 9 < 2 the function N x ix) is monotonously increasing. 
In the stationary case the presence of a decreasing branch 
implies that different values of x ar e associated to the 
same Cf/C|. 

A stationary solution of the diffusion equation Eq. I|17|l . 
fulfilling the boundary conditions Eqs. I|19l) . reads 



a D l 
giving the surface concentration 



Cf 



Gmix) 



(26) 



(27) 
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FIG. 3: Stationary functions N x (x) and Na(x) vs - the solid 
state composition x- The chosen parameters are 9 = 3, 
f^e^ = 1, and Gi = Gi. For this choice the stationary value 
X st is 0.5. The dashed curve illustrates the typical construc- 
tion of the limit cycle for the developed instability. 



The simplest form of a possible time evolution giving 
rise to oscillatory behavior is shown in Fig. [31 If there 
is a time-scale separation between the x-variations and 
the variations of the surface concentrations Cf the system 
would move from the left maximum of N x (x) to the right 
until x ~ (^A^XCf/Cf). Then the system has time to 
adapt (Cf/Cf) to the present value of x thereby moving 
down along the A^ x (x)-curve. Once the the minimum 
is reached (d/dt)(C( /Cf) will drive the trajectory away 
from N x (x) quickly resulting in |x| > \(d/dt){C{ /Cf)| 
which will move it towards the left branch of N x (x). 
Then the second half of the oscillation may start. As 
will be shown below via numerical simulations the ac- 
tual behavior is somewhat different. In any event, the 
general possibility of oscillatory behavior requires a non- 
monotonous behavior of N x (x), i.e. 9 > 9 C = 2. How- 
ever, the linear stability analysis will reveal that this con- 
dition is not sufficient. 



B. Linear stability analysis 



1) 



This yields (using again /^e' 

in the stationary limit. The function Nq(x) is also shown 
in Fig. |21 Both the functions N x (x) and Nq (x) intersect 
at x — X st given by 



X 



G\ + G2 



Gi 
G 



(29) 



thereby defining the total incoming flux G — G\ + G 2 - 
It defines the stationary point of the present system for 
which both Eq. (|17l) and Eq. (|21|l have vanishing time- 
derivatives. 

Now we discuss some immediate consequences of the 
properties of the functions N x (x) and Ng(x)- For this 
purpose we consider the (x, Cf/Cf) -plane. These func- 
tions separate different regions of the phase plane. Natu- 
rally, a point on this plane does not describe the complete 
system configuration because the diffusion is reduced to 
C(z = 0). For Cf/Cf > N x (x) the time derivative of 
x(t) is positive and vice versa. This means that a system 
above N x (x) is driven toward larger \, whereas a system 
below behaves oppositely. 

In case of a standard relaxation oscillator the equa- 
tion for x would be complemented by an equation for 
(d/dt)(C( /Cf). Such an equation does not exist for the 
present model. However, an implicit time evolution of 
(Cf/Cf) is expressed by Eq. (J2|. When (d/dt)x(t) is 
very small the diffusion field can adjust to the boundary 
conditions. This way the ratio Cf/Cf will be adjusted 
until Eq. J2H1) is fulfilled. In particular, Cf/Cf will de- 
crease for (Cf/Cf) > Ng(x) an d vice versa. This is 
reflected by the arrows in Fig|3 



The stability of the stationary crystal growth is an- 
alyzed with respect to infinitesimal perturbations of the 
solute distribution in the solution and the solid state com- 
position around x = X st and the corresponding values of 
Cf, given by Eq. which will be denoted Cf. We 
choose 

6Ci(z, t) = 8Ct exp (7* - Pi z) (30a) 

and 



$x(t) =Sxexp (jt). 



(30b) 



Here 8Cf and $x are the amplitudes whereas the com- 
plex wave number pi = Kepi + i Irapi describes the decay 
of the concentration perturbations above the crystal sur- 
face. It requires Kepi > 0. The instability arises when 
the real part of the perturbation increment 7 becomes 
positive. The chosen form (|30|l is compatible with the 
time-evolution close to the stationary point. 

The infinitesimal form of Eq. (|21|l is of primary inter- 
est. A short calculation gives 



^ = -i + x st (i-x st ) 



26 



SCi 



5C 2 



CfS x GfS x 



(31) 



The relevant ingredient 6Cf / (CfSx) can be obtained 
from the boundary condition l|19fl as 



8Cf 



ad In Ti I d\ 



Cf5 X a + D iPl T t {x s 
Finally, from Eq. (|17|l one has 

7 = Dip\ = D 2 pl , 



(32) 



(33) 
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i.e. the relation between the pi and 7. Combining 
Eqs. i|5T)l. ijSDl and and using the specific depen- 
dencies (|20|) we finally get 



-2^ = -i + x st (i-x st ) 

a £ G 



{0 + 4>) 



D1T1P1 + 



] D2T2P2 
D2T2P2 + 



(34) 



How does the nature of the stability depend on the 
total flux G? One always has Cf oc G. Furthermore, in 
the limit of large G one also obtains 7 oc G and, using 
Eq. pi oc \fG. As a consequence one has SCi oc 

Ci/pi oc VG. Thus, Eq. IpTTJ) boils down to 



7 = a 2 G[-l + X st (l-X st )20] 



(35) 



C. Exact solution 

From Eq. I|34|l the critical value G c can be determined. 
As shown in the Appendix the critical value G c in para- 
metric form is given by the expression 



G c 



1 



(i-2 X )-e 



(0 + 0)A (0-0)/A 



(36) 



.(V2CA + 1) 2 + 1 (V2C/A + 1) 2 + 1 
Here the variable £ is the root of the equation 

(fl + <^(V2CA)+(fl-^(V2C/A) = x{l \ x) , (37) 

where the function is defined as 

x(x + 1) 



V ; (a; + 1) 2 + 1 
and the function A(x) is 

A(x) = A Q . 5 e e ^» 2 
with the prefactor A0.5 given by the expression 




(38) 



(39) 



(40) 



Since <j> < 9, both terms on the left-hand side of 
Eq. 1|37[) are non-negative increasing functions of the ar- 
gument £. The maximum of their sum is equal to 29 
whereas the minimum of the right-hand side at % = 0.5 is 
equal to 4. So, we obtain again the condition 9 > 9 C nec- 
essary for instable behavior. In this case for 2#x(l — \) > 
1 equation l|37|l possesses only one root which together 
with |nni> dete rmine the lower boundary G c of the in- 
stability region. In the close vicinity of the threshold, 
< 2#x(l - x) - K 1 (for x ~ 0.5) this equation can 
be solved analytically, giving the expression 



{(9 + 0)/A O r o + (9 



)A , 



4V2[20 X (1 - x ) - 1] 



(41) 



One has 7 > exactly when x st € (x-,X+), i-e. X s * is and thus 
on the unstable branch of N x (x)- In contrast, for G — > 
one has 7 = —a 2 G < 0. This can be seen from Eq. 134|) 
because in this limit also pi — > 0. Thus there exists a 
critical flux G c such that Re 7 = for G = G c and Re 7 > 
for G > G c - Since 7 = pi = cannot be a solution 
from Eq. (|34J) the disappearance of Re 7 implies that 7 is 
purely imaginary. 

Our goal is, first, to determine G c explicitly and, sec- 
ond, to understand its dependence on the model pa- 
rameters on a more qualitative level, i.e. G c — G c (x I 



[(9 + ^)/A . 5 + (9 



)A . 5 ] 2 



[20*(1 - X) - 1] £ 



(42) 



For the system with <f> = the asymmetry of the sol- 
ubilities substantially increases the critical diffusion flux 
G c . Indeed, since the species diffusivities in solutions are 
typically of the same order, D\ ~ D2, the species sol- 
ubility difference reflected in the coefficient f v -^i 1 or 
/, ; ^> 1 matches (see Eq.l^UJl) the inequality A .5 ^> 1 or 
A0.5 -C 1, respectively. For « 9 the term (9 — 0)Ao.5 in 
(|42|l can be neglected and an increase of A0.5 gives rise 
to a remarkable decrease in G c - 

Equation l|37|l was solved numerically to analyze the 
system behavior far from the threshold 9 C = 2. For 9 = 3 
the results for the x-dependence of G c are presented in 
Fig.0] In addition we have included N x (x) in Fig.^] The 
upper frame exhibits the results for = 0, whereas the 
lower one contains the strongly asymmetrical case with 
(f> = 9 = 3. 

The symmetrical system exhibits minimal G c when the 
species have the same solubility. A difference in solubility 
as reflected by the change of A0.5 by a factor of ten causes 
G c to increase by a similar factor. For such A0.5 the 
G c (x)-curves become asymmetrical with respect to x — 
0.5. Exactly this case should be characterized by the 
instability forming the limit cycle constructed in Fig. [31 
following the standard notions of relaxation oscillations. 

The system behavior for <j) ~ 9 is distinctly different as 
can be seen in the lower fragment of Fig.0] In particular, 
for 4> = 3 a ten-fold increase in A0.5 induces a hundred- 
fold drop in G c . It should be noted that A0.5 = 1 corre- 
sponds to 77 = <p/2 which is non-zero here. Furthermore, 
for A0.5 = 10 the dependence G c (x) is highly asymmet- 
rical and passes many orders of magnitude. For large 
values of A0.5, the left part of the decreasing branch of 
N x {x) can be unstable whereas the right half can be sta- 
ble. In this case the limit cycle of the developed oscilla- 
tions deviates substantially from the classical form. 
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FIG. 4: The critical value of G c vs. x f° r # = 3 and some 
different A0.5 representing the difference in the species solu- 
bilities. The upper frame corresponds to a symmetric system 
with <f> = 0, the lower one to a strongly asymmetrical one with 
cj> — 3. The dashed curves correspond to the function N x (x)- 



interesting to analyze the impact of the asymmetry <j>. In 
what follows we will restrict ourselves to the case <fi = 6, 
where SC2 = and the damping is due to the concentra- 
tion fluctuations of species 1. Changing from <fr = to 
4> = 9 the relevant term in Eq. I|31|l . characterizing the 
concentration fluctuations and thus the reduction of the 
real part of 7 changes from — (I/T1+I/T2) to — 2/t\ which 
corresponds to a change by a factor of 2/(1+ti/t 2 ). Note 
that 



r 1 /r 2 = A2. 5 exp[0(l-2 X st )] 



(45) 



When A0.5 ^ 1 the result of the concentration fluctu- 
ations is thus reduced which has a positive cumulative 
effect on the instability onset and, as a result, G c should 
be decreased. Moreover, this ratio increases as x de- 
creases. This rationalizes the asymmetry of the G c (in- 
dependence with smaller values located on the left-hand 
side. Exactly these features emerge from the solution of 
the exact equation 1|34|> . 



V. NUMERICAL ANALYSIS OF THE 
NONLINEAR DYNAMICS 

When the growth becomes unstable nonstationary pat- 
terns in the solution and the induced pattern in the crys- 
tal bulk develop. In order to analyze their characteristic 
properties the system of equations in Sect. ITTT1 was stud- 
ied numerically. First, we introduced the spatial and 
temporal scales l sc and r sc 



\ a 2 



(46) 



and, in addition, a parameter having the dimension of 
concentration 



D. Instability mechanism: qualitative description 

It is possible to obtain a better understanding of how 
the degree of instability depends on the system parame- 
ters. We use sufficiently large G such that a <C \DiPin\. 
The values of t% are always analyzed at x = X st - F° r 
reasons of simplicity we also assume D\ = D2 = D (im- 
plying p = pi = P2)- Then we can rewrite Eq. ll.'il'l) as 



SCf 



a dhiTi/ d\ 



CfS X pD 



thereby 
SCf 



5C S , 



Cfdx CfS X pD 



Tl 



T2 



Tl 



Tl 



(43) 



(44) 



The cumulative effect of these terms gives rise to a de- 
crease of the real part of 7. Thus, concentration fluctua- 
tions always tend to stabilize the stationary point. It is 



2T 



(47) 



Then we rescaled time t and spatial coordinate z as well 
as the species concentrations Ci in these units 



t 



t, 



Ci 



G sc • Ci 



For the sake of simplicity we have kept the same desig- 
nations for these variables. The fluxes Gi, G were also 
measured in units of 

1 

Gar, = 



VD^D~ 2 t 2 ' 



namely, 



{G, Gi) — > G S c • (G, Gi 



(48) 



Then the obtained system of dimensionless equations was 
integrated using the Crank- Nicholson method. The tem- 
poral and spatial steps of integration as well as the system 
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size L were chosen such that the dynamics be practically 
independent of their particular values. 

The following four specific cases were analyzed to 
demonstrate characteristic features of the physical sys- 
tem. The first two ones are the symmetrical model with 
parameters 9 — 3, <fr — 0, A .5 = 1 and Xst = 0.5 at 
the initial stage. They differ in the total diffusion flux 
G = 1 and G = 50. The case with G — I describes the 
system dynamics not too far from the instability bound- 
ary G c « 0.6 (see Fig. QJ. Under these conditions the 
instability was expected to demonstrate quasiharmonic 
behavior. The case G = 50 corresponds to the substan- 
tially nonlinear stage of the instability. 

The other pair of cases are the asymmetrical model 
with = 3 and = 3 being actually a limit situation that 
can be considered accurately within the present analysis. 
To single out the effects caused by nonzero values of cj) 
we restrict ourselves to A0.5 = 1 and Xst = 0.3. Cor- 
respondingly, G c {x) is minimal for this value of x- The 
values of G are set to 0.5 and 5. 

Figure [5] visualizes the surface dynamics for the first 
two cases. As seen in the frames of its left column, the 
oscillations of Cf and x are rather harmonic as expected. 
The resulting limit cycle is shown in the lower left frame. 
Again, as expected, this cycle is of elliptical form located 
along the unstable branch of the nullcline N x (x)- Even 
for G = 1 which is close to G c ss 0.6 the oscillations are 
already determined by the nullcline. 

For G = 50 the dynamics become relaxation-like and 
X plays the role of the fast variable. As seen in the right 
column of Fig. [5] the time pattern x(i) consists of a se- 
quence of slow motion fragments joined by rather sharp 
jumps. In agreement with classical relaxation oscillations 
the fragments of slow dynamics correspond to the sys- 
tem motion along stable branches of N x (x) whereas the 
sharp jumps describe the fast transitions between these 
branches. However, the surface species concentrations 
exhibit anomalous behavior from the standpoint of clas- 
sic relaxation oscillations. Most importantly, Cf/C| dis- 
play some remarkable variations during the fast dynam- 
ics phase of x- This is still reminiscent of the result of 
the linear stability analysis where SCf cx 6x- As a re- 
sult the lines of the limit cycle connecting the increasing 
branches of the nullcline N x (x) are not horizontal but 
inclined to the x _ axis at a certain visible angle. It should 
be noted that this inclination is not a standard conse- 
quence of the finite ratio between the time scales of fast 
and slow dynamics because, otherwise, the form of these 
curves would deviate form the straight lines substantially. 
It is due to the existence of the simultaneous change in 
both the species concentrations during the period of fast 
motion. 

The asymmetrical system with 9 = 3, <fi — 3 is depicted 
in Fig. HI] for G — 0.5 and G = 5. The basic feature dis- 
tinguishing this system from the symmetrical one is the 
strong asymmetry of the critical diffusion flux with re- 
spect to x — 0.5. In particular, as shown in Fig. ^] for 
G = 0.5 only the left half of the decreasing branch of the 



nullcline N x (x) is unstable whereas for G = 5 the insta- 
bility region is located approximately within the bound- 
aries 0.2 < x < 0.6. The resulting limit cycles are shown 
in the lower frames of Fig. EI F° r G = 0.5 the stationary 
point is not too far from the instability boundary and the 
limit cycle is located in a rather narrow neighborhood of 
the unstable branch of N x (x)- For G = 5 the separation 
of the system dynamics into the slow and fast motion 
fragments becomes pronounced. In this case the limit 
cycle embraces even an increasing branch of the N x {x)- 
However, its lower part after passing the minimum near 
X = 0.8 continues to follow the formerly unstable frag- 
ment of the decreasing nullcline branch until it reaches 
the instability boundary at x ~ 0.6. Only after passing 
this point it leaves the branch and jumps to the oppo- 
site stable branch of N x (x)- This effect of "adhesion" to 
the unstable branch of N x (x) is of another nature than 
the well known "French duck" fragment of the limit cycle 
for standard relaxation oscillations (see, e.g., Ref. |27j1. 
caused by the close proximity of the stationary point to 
the extreme of the nullcline. 

The two time patterns of the observed oscillations de- 
viate substantially in shape both from the quasiharmonic 
oscillations and relaxation oscillations. Especially for the 
flux G = 5 the found pattern looks like a sequence of pro- 
nounced spikes joined by fragments of slow motion along 
the left decreasing branch of N x (x)- 

The resulting crystal profiles are shown in Fig. \7\ Ba- 
sically, the profiles are scaled mirror images of the time- 
dependent x as shown in Figs. [5] and EJ because the 
growth velocity only weakly varies with time. As should 
be expected the spatial oscillation period is orders or 
magnitude larger than the atomic scale. 



VI. DISCUSSION 

We have presented a ID model for OZ. The growth 
rate as the central non-linear coupling term between the 
solution and the crystal was derived based on layer-by- 
layer crystal growth mechanisms. For this purpose parti- 
cle adsorption, surface diffusion and finally desorption or 
inclusion processes at the steps are taken into account. 
This way, adsorbed particles do not only experience the 
local environment of their adsorption site, but a much 
more averaged one depending strongly on the composi- 
tion of the crystal surface. This is an essential feature 
of the present model, because together with volumetric 
species diffusion it may provide the synchronization ef- 
fect necessary to successfully describe the experimental 
findings in more than one dimension. 

Our model (Eqs.JT7J|-|J2IJ) differs in several respects 
from the existing model of L'Heureux et al. (i) 
The most essential difference is the choice of the growth 
mechanism, i.e. the definition of the growth rate (%) . 
In Refs. [f5|, [l6( a phenomenological equation based on 
local atom adsorption is formulated. It is one of the 
simplest polynomial expression for fi(x) which displays 
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FIG. 5: Dynamics of the symmetrical system 8 = 3, (f> = 0, A0.5 = 1, Xst = 0.5. The left column has the results obtained for 
G — 1 which is larger than G c = 0.6 , see Fig. 2] The right column exhibits the results for G = 50. 



nonlinear behavior, (ii) The boundary condition at z = 
used in |15l Il6| for crystal growth from solution can be 
written as 

D i dC i /dz = Xi{n + r 2 ) (49) 

in our notation. Here one uses xi = 1 — X2 — X- It is 
identical to Eq. I|19|) if (d/dt)x(t) = 0. But as soon as 
X changes with time both boundary conditions are not 



identical. Actually, to derive Eq. 14911 one can formu- 
late the mass balance in a newly generated crystal layer. 
Then Eq. I|49() only emerges if the x~ va ^ ue i n the layer at 
some time, determining the nature of the growth rates 
ri{x) during the growth of the new layer is identical to 
the %-value in the resulting new layer. In general this is 
not the case. Thus we feel that Eq. (|19fl is more gen- 
erally valid although from a practical point of view no 
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FIG. 6: Dynamics for the asymmetrical system 6 — 3, (f> = 3, A0.5 = 1, Xst = 0.3. The left column presents the data obtained 
for the total diffusion flux G = 0.5, whereas the right one exhibits the same data for G = 5. 



essential difference should be present, (iii) We specify 
the gradient at z = L whereas in [TEl Ht| the concen- 
tration is given. Our choice is conceptually more simple 
because the value of L does not enter the calculations (see 
e.g. Eq. (|26p). However, the underlying physical picture 
does not change whether the gradient or the concentra- 
tion itself are fixed as long as L is sufficiently large, (iv) 
Another distinctive difference is the explicit calculation 



of the diffusion field in the solution above the crystal. We 
do not apply the boundary layer approximation, regard- 
ing a crystal growing through a supersaturated solution 
[l5j, but we neglect the spatial growth with respect to 
the diffusion. This way, we obtain a source-sink system 
with a diffusion field that has to be treated explicitly, 
because the interaction of particle accumulation or de- 
pletion on the solution side and the autocatalytic growth 
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FIG. 7: Spatial patterns formed in the growing crystal. 



on the crystal side is essential for the existence of OZ. In 
addition, we perform the stability analysis without fur- 
ther approximations. This may be essential because the 
oscillatory dynamics can stem from the counteraction of 
strong "forces" mutually compensating each other at the 
first approximation. We would like to mention that in 
|16| the system of equations has been numerically solved 
without invoking the boundary layer approximation. On 
a qualitative level similar results were obtained as com- 
pared to the analytical treatment within the boundary 
layer approximation, (v) The surface roughness param- 
eter, used in |l5j| . is not necessary in the present analysis. 

The proposed growth mechanism is valid for a <C l s <C 
I. In the limit of l s ~ a incorporation would be governed 
by the local crystal composition, again. In case of l s « I 
desorption can be neglected and nearly every adatom will 
be incorporated regardless of its type and \. This would 
result in crystals exhibiting the stoichiometric composi- 
tions |2^| of the influx. 

The present model cannot exhibit the bistability found 
in [l^, because the composition of the only stationary 
point is determined by the influx ratio. Any crystal 
growth with a composition different from this ratio will 



result in a buildup of the currently "disfavored" compo- 
nent until its growth rate increases. Thus the oscillations 
have to revolve around or run into this fixed point. This 
corresponds well to the model picture, where a bistabil- 
ity is only possible if either the supersaturation is not 
high enough with respect to the minor component or if 
complementary crystals grow in close vicinity. 

The time scales of diffusion and crystal growth define 
the qualitative behavior of the oscillations. At very low 
concentrations the crystal growth rate and consequently 
the changes in crystal composition are very slow. The 
diffusive processes on the other hand are fast enough to 
counteract this and no oscillatory behavior is observed. 
With rising solute concentrations the growth rate in- 
creases as well, whereas the characteristic speed of dif- 
fusion stays constant resulting in soft transitions and si- 
nusoidal oscillations. At even higher concentrations, the 
speed of crystal composition change supercedes the dif- 
fuse reaction of the solution by many orders, creating 
sharp transitions, closely following N x (x)- This qualita- 
tive description coincides with the results from the linear 
stability analysis At low influx it is stable possessing two 
imaginary eigenvalues with negative real parts. At G c the 
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real part changes sign and the system becomes unstable. 

Finally, the numerical calculation of asymmetrical sys- 
tems show the stabilization of a formerly unstable branch 
of N x (x) ■ If the fixed point is shifted far enough onto the 
stable branch, the trajectories slowly approach it along 
the nullcline, but upon nearly reaching it are directed 
away with the onset of new oscillations. The amplitude 
of these oscillations can have continuously growing char- 
acter or they can exhibit certain characteristic values. 
This might be the onset of frequency doubling and pro- 
vide a route into chaotic behavior of the system. This 
however, is beyond the scope of the present manuscript. 



VII. CONCLUSION 

The presented boundary reaction diffusion model, de- 
scribing OZ of crystals grown from solution. Thereby it 
applies a mechanism including surface diffusion, which 
can be derived from microscopic properties. OZ arises 
due to diffusive build up of the disfavored component 
above the crystal until the autocatalytic growth process 
in combination with the large population turn the crystal 
composition in favor of the formerly disfavored species. A 
linear stability analysis was carried out without further 
approximations. The results correspond to the experi- 
mental findings of a supersaturation threshold before the 
onset of OZ. Numerical simulations of asymmetric cases 
show a stabilization of previously unstable parts of the 
system. 
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EIGENVALUE EQUATION FOR THE 
BOUNDARY-REACTION DIFFUSION MODEL 

Let us convert the system of equations and iTHI) 
into a form more appropriate for its analysis. By virtue of 
(|33fl both the wave numbers pi have the same argument 
ip € (— 7r/2; it/2), i.e. pi = \pi\ exp(iip) whose twice value 
gives the argument of 7 and their absolute values obey 

1 /2 

the equality \p t \ = (M/A) ■ 

A new variable £ > and parameter A introduced as 
follows 



V-D1-D2T 2 _0(i_ 2x )+0| 




4>+8(l-2 x ) 

/ D 2 r 2 V D2' 
enable us, first, to write 

Dm \ Pl I = <A , D 2 t 2 \ P2 1 = </ A 
and, then, to represent equation l|34|l in the form 

(C/A)e»* 



(50) 
(51) 

(52) 



-i + x(i-x) 

(CA)e^ 
' (CA)e i " 



1 



{e~<p) 



(C/A)eW + 1 



(53) 



where the parameter g stands for the dimensionlcss 
species flux casing the crystal growth 
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^DjT 2Tl T 2 G = VA-DarV^-^+^G . (54) 



Finally the split of equality (|53|l into the real and imagi- 
nary parts yields the desired coupled equations specifying 
actually the eigenvalue 7 



2-cos(2^) = -l + x(l-x) 
9 

c 2 

— sin(2-0) = x(l — x) sin ip 
9 



(CA)(CA + cos^) 
(CA + cos?/') 2 +sin 2 ^ 

m ^. 

(£A -I- cos ip) 2 + sin ip 



(C/A)(C/A + cos?/*) 
(C/A + cos-0) 2 +sin 2 '0 
, (C/A) 
(C/A + cos?/0 2 + sin 2 i/;_ 



(55a) 
(55b) 



As follows from equations l|55|l the value ( = is not a 
root of this system. So the considered instability is to 
arise through the real part of 7 changing the sign with 
its imaginary part having some finite value. Therefore 



to find the instability boundary Re 7 = we can set 
ip = 7r/4 in the given equations. In this way the former 
equation l|55a|l converts into one specifying the value of 
C for chosen values of the parameters x? $j 4>i an d A 
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(V2CA)(V2CA + 1) 



.(>/2C/A)(V5C/A + l) 



v TJ (\/2CA + 1) 2 + 1 v r; (V2C/A + 1) 2 + 1 X(l-X)' 
whereas the latter equation determines the critical value of the diffusion flux g c 



(56) 



V2\(l - X) 



+ . 



1/A 



(V2CA + 1)2 + 1 



(V2C/A + 1)2 + 1 



(57) 



Returning to the dimensional variables we get expressions l|36l) and l| 37|). 
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